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Gravitational perturbations and metric reconstruction: Method of extended 
homogeneous solutions applied to eccentric orbits on a Schwarzschild black hole 

Seth Hoppc iE and Charles R. EvansQ 
Department of Physics and Astronomy, University of North Carolina, Chapel Hill, North Carolina 27599 

We calculate the gravitational perturbations produced by a small mass in eccentric orbit about 
a much more massive Schwarzschild black hole and use the numerically computed perturbations 
to solve for the metric. The calculations are initially made in the frequency domain and pro- 
vide Fourier-harmonic modes for the gauge-invariant master functions that satisfy inhomogeneous 
versions of the Regge- Wheeler and Zerilli equations. These gravitational master equations have 
specific singular sources containing both delta function and derivative-of-delta function terms. We 
demonstrate in this paper successful application of the method of extended homogeneous solutions, 
developed recently by Barack, Ori, and Sago, to handle source terms of this type. The method 
allows transformation back to the time domain, with exponential convergence of the partial mode 
sums that represent the field. This rapid convergence holds even in the region of r traversed by the 
. point mass and includes the time-dependent location of the point mass itself. We present numerical 

results of mode calculations for certain orbital parameters, including highly accurate energy and 
angular momentum fluxes at infinity and at the black hole event horizon. We then address the issue 
of reconstructing the metric perturbation amplitudes from the master functions, the latter being 
weak solutions of a particular form to the wave equations. The spherical harmonic amplitudes that 
represent the metric in Regge- Wheeler gauge can themselves be viewed as weak solutions. They are 
in general a combination of (f) two differentiable solutions that adjoin at the instantaneous location 
of the point mass (a result that has order of continuity C _1 typically) and (2) (in some cases) a 
delta function distribution term with a computable time-dependent amplitude. 

o : 

O" 1 PACS numbers: 04.25.dg, 04.30.-w, 04.25. Nx, 04.30. Db 

u . 

(N ' I. INTRODUCTION 

> 

Considerable research on the two-body problem in general relativity has been fostered over the past decade by the 
Q^ • prospects of detecting gravitational radiation from extreme-mass-ratio binaries. The general relativistic two-body 
problem is notoriously difficult, as it involves dynamics of the motion of the bodies and of the gravitational field 
\q , itself. Gravitational wave emission carries away energy and angular momentum from the orbit, leading to inspiral and 
eventual merger. The future joint NASA-ESA LISA mission [l| is expected to detect between tens and thousands of 
such extreme-mass-ratio inspirals (EMRIs)-binaries composed of a compact object (fj, ~ 1 — 50M Q ) in orbit about a 
~~ * 1 supermassive Kerr black hole (M ~ 10 5 — 10 7 M Q ) out to cosmological distances (z ~ 1) |2J. The small mass ratio 
10~ 7 < fi/M < 1CP 3 of expected astrophysical sources [3[ implies a gradual change in orbital parameters, with > 10 5 
• '-J i wave periods as the binary evolves through the LISA passband (10~ 4 — 1CP 2 Hz). Detailed theoretical calculations will 
aid in both detection of EMRI gravitational wave signals and in determination of the source's physical parameters. 

Quite apart from the prospects of astrophysical observation, this problem is one of intrinsic interest in theoretical 
physics. Of the various possibilities, the physically simplest compact binary is one composed of two black holes. 
Such a system eliminates the complications of stellar microphysics and reduces the problem to a minimum parameter 
set. In approaching the problem mathematically, the extreme mass-ratio and gradual orbital evolution is of benefit 
theoretically, allowing black hole perturbation theory to be used. Furthermore, the small mass ratio allows even the 
black hole structure of the small mass to be ignored (at lowest order), restoring a point-like (particle) behavior [4j on 
length scales that are large compared to /x and thereby simplifying the perturbation problem. 

The perturbation problem proceeds in stages. At the outset the motion of the particle is taken as a geodesic 
(/i/Af — > 0, or zeroth order) on the background spacetime. The first-order (in n/M) gravitational field perturbation 
is then computed, yielding a new metric g^ v = g^ v +Pfj,u that corrects the background metric g^. The gravitational 
waves in the perturbation p^ v carry energy and angular momentum to infinity and down the black hole event horizon, 
giving rise to a back reaction or local self-force (SF) on the particle that has both conservative and dissipative terms. 
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Formally, the SF depends on gradients oip^v and acts locally on the particle to accelerate it off its background geodesic. 
Once the first-order correction to the motion is successfully computed, the calculation may proceed to second order 
in the field perturbation (see Pound [5] for a recent background discussion and an alternative formulation). 

Yet having idealized the small body as a point particle, the metric perturbation and SF are found to diverge at 
the location of the particle, and the formal perturbation to the equation of motion is meaningless without careful 
regularization. This problem is similar to the classic SF problem of an accelerating, radiating charge in electromagnetic 
theory in flat spacetime [f|. Two pivotal papers, by Mino, Sasaki, and Tanaka p} and Quinn and Wald @, showed 
how the metric perturbation may be separated into a divergent, direct part p^" and a finite tail term p^£r, with the 
latter providing the regularized field that makes the SF finite. As an alternative, Detweiler and Whiting Q proposed 
decomposing the metric perturbation into regular pff v and singular p^ u parts. Under this interpretation, p& is a 

solution to the vacuum field equations, but gives rise to the same SF as P^f - 

Since then, SF calculations have been made in certain special cases [i0l - [T3 |. See the review by Barack Q. Ul- 
timately, the theory aims to provide self-consistent SF calculations of arbitrary orbits about Kerr black holes. In 
this paper, we concern ourselves with a more modest goal: demonstrating a complete computation of the radiative 
gravitational perturbations produced by a mass in eccentric orbit on a Schwarzschild black hole and reconstruction 
of the corresponding parts of the perturbed metric in Regge- Wheeler gauge. While we leave for another occasion 
computation of both the nonradiative perturbations and the SF, the accurate reconstruction of the radiative parts of 
the metric, at all locations up to and including the point mass, should serve as a starting point for a further gauge 
transformation or alternative regularization technique. 

We note in passing that most work to date computing EMRI evolution has not made use of local SF calculation. 
Sufficiently adiabatic changes in an orbit on Schwarzschild spacetime allow a balance calculation approach [l5| , where 
orbital energy and angular momentum are "evolved" (acausally) to match corresponding gravitational wave fluxes 
through bounding surfaces at large radius and near the horizon. Much effort is ongoing to extend the reach of 
adiabatic calculations [l6l - [l8j . Unfortunately, the approach only approximates dissipative SF terms and cannot 
account for conservative SF effects. In any event, the more self-consistent SF approach should serve to confirm the 
validity of these or other approximations. 

Perturbation theory for Schwarzschild black holes has a traditional formalism pioneered by Regge and Wheeler [l9| , 
Zerilli [20], and Vishveshwara 21 1 that uses spherical harmonics and the Regge- Wheeler gauge to simplify algebraically 
the form of the metric perturbation. At each spherical harmonic order there are just two master functions, \&^ n (t, r) 
and ^1^(t,r), one for each parity or gravitational degree of freedom, which satisfy linear inhomogeneous wave 
equations in t and r. The formalism was improved by Moncrief (22| and colleagues [23|, [24|, makinguse instead of 
gauge-invariant master functions that satisfy similar wave equations. Recently Martel and Poisson [25[ have placed 
the theory in both a gauge-invariant and covariant form. 

For perturbations of Kerr black holes, Teukolsky [26[ developed a formalism based on Newman- Penrose curvature 
scalars and spin-weighted spheroidal harmonics. In the frequency domain the radial part is a single (complex) master 
equation [27| , which can, of course, be applied to a Schwarzschild hole as well [HI, H|| . 

An alternative to the Regge- Wheeler- Zerilli (RWZ) approach has recently been advanced by Barack and Lousto [29| . 
They propose directly evolving the ten spherical harmonic amplitudes that describe the metric perturbation in Lorenz 
(or harmonic) gauge. In this direct metric perturbation approach, the equations separate into even- and odd-parity 
sectors, yet still involve systems of seven and three coupled equations, respectively. Barack and Sago [ill, have 
used the formalism to compute the time evolution of metric perturbations generated by circular and eccentric orbits 
on Schwarzschild, along with the resulting SF components. 

The RWZ and direct metric perturbation approaches each have advantages and disadvantages. The direct metric 
perturbation formalism yields directly what one wants as an input to a SF calculation, namely the metric itself in 
Lorenz gauge. In a time domain calculation, as so far employed, it has the disadvantage of requiring simultaneous 
solution of a large set of coupled partial differential equations (PDE's). Anticipating the subtraction involved in the 
SF regularization, Barack, Lousto, and Sago have built a fourth-order convergent finite difference code to compute 
the modes to sufficient accuracy. In contrast, the RWZ approach has the advantage that only a single uncoupled 
wave equation need be solved for each mode and parity. Unfortunately, an added step is required to reconstruct the 
metric from the mode solutions. Moreover, the reconstruction involves terms that are singular at the particle location 
and the simplest reconstruction yields the metric perturbation in Regge- Wheeler gauge [3(| HH . Finally, the RWZ 
approach provides only the radiative (£ > 2) parts of the perturbation and the nonradiative modes (£ = 0, 1) must be 
derived by separate means. 

In this paper we opt for using the gauge- invariant RWZ approach detailed by Martel and Poisson [25|, and adopt 
the Zcrilli-Moncrief ^™ = ^fm" and Cunningham-Price-Moncrief X I'^' M = master functions for even and 

odd-parity, respectively. Our use of this relatively standard method is augmented, though, by a new technique that 
enables accurate reconstruction of the corresponding parts of the metric in Regge- Wheeler gauge. We leave for a 
later occasion our own consideration of the monopole and dipole terms (which are essential to a SF calculation) and 
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instead direct attention to discussion by Detweiler and Poisson [32| and recent successful numerical implementation 
by Barack and Sago [l4j . 

The master functions can be obtained directly by numerical evolution (solution of PDE's) in the time domain (TD) 
(see e.g., [13, HH, EH EH US l33T - l35j ]) or by numerical inte grat ion of ordinary differential equations (ODE's) for the 
Fourier modes in the frequency domain (FD) (see e.g., (la. I36l438l |). Each method has strengths and weaknesses. TD 
calculations require solving just one equation for each £, m mode and time dependence of the subsequently recon- 
structed metric and SF is of direct interest. Disadva ntag es of TD calculations include (1) modeling the discontinuous 
source movement through the finite difference grid [1J, l29[; (2) numerical stability of PDE evolution; (3) difficulty 
devising numerical schemes of adequately small truncation error; and (4) challenges in posing outgoing wave boundary 
conditions at finite radius. In contrast, in FD calculations (1) the numerical errors tend to be much smaller (i.e., by 
solving ODE's); (2) outgoing wave boundary conditions are handled mode-by-mode and extrapolated to infinity and 
to the black hole event horizon; and (3) the discontinuous source presents few difficulties in computing (at least) the 
Fourier mode functions Re m n(r)- However, FD methods require, for eccentric orbits, computing and summing over 
numerous harmonics n of the radial libration frequency Q r for each £, m and transformation to the TD is nontrivial 
given the singular source terms. 

Barack, Ori, and Sago (BOS) [38] highlighted the latter difficulty. They used the model problem of a scalar field 
$(£, r, 9, if) generated by a scalar point charge in eccentric orbit on Schwarzschild. The spherical harmonic modes 
<j>e m (t,r) = r$e m (t,r) satisfy a wave equation with a singular source, Sf^ iaT (t, r) — Cg m (t,r) S[r — r p (t)]. Here 
Ce m (t,r) is some smooth function and r = r p {t) describes the radial libration of the particle's worldline between two 
turning points. In the FD, ODE's are solved for the Fourier-harmonic modes Re m n(r). These mode functions are, 
at each point r, Fourier series coefficients. The resulting Fourier series converges for the piecewise continuous (C°) 
4>im(t, r) but the singular nature of the source S makes </W(t, r) converge slowly in the region traversed by the point 
charge. The radial derivative d r ^>e m is however discontinuous at r — r p (t) and its Fourier series only converges, in the 
usual sense [39j . almost everywhere. The attempt to assemble the radial derivative from the Fourier series is plagued 
by the Gibbs phenomenon; the series converges to the mean value at the discontinuity and the series "overshoots" 
and fails to converge properly in the limit as both n — > oo and r — > r p (t)^ . 

BOS circumvented the difficulty with a new method of extended homogeneous solutions. In brief, they use FD 
analysis to find Fourier-harmonic mode solutions to the homogeneous equation, valid outside and on either side of the 
source libration region. The associated Fourier series converge exponentially fast to homogeneous solutions of the TD 
wave equation. They then analytically extend both homogeneous TD solutions into the source libration region up to 
the instantaneous position of the point charge. Summed to adequately high order, the two homogeneous solutions 
match in value at r p (i), as expected. With the field represented in this way, the left and right derivatives can be 
accurately determined. BOS argued that the method should work for other problems with similar wave equations, 
including the Teukolsky equation. 

We show in this paper that the method can indeed be extended to the case of gravitational perturbations computed in 
the RWZ formalism, and apply the method to a large set of Fourier- harmonic modes stemming from a mass in eccentric 
orbit on Schwarzschild. (Note that Barack and Sago fl4j previously implemented this method in the gravitational 
case but only for the £ = 0, 1 modes in Lorenz gauge.) An important distinction arises: in the gravitational case the 
source distribution in the Regge- Wheeler gauge contains both delta function and derivative-of-delta function terms, 

S im (t, r) = G lm {t, r)5[r~ r p (t)} + F im {t, r) S'[r - r p {t)}, (1.1) 

with Ge m (t, r) and Fi m (t, r) being smooth functions. As a consequence the master functions have a jump discontinuity 
at r = r p (t) (referred to sometimes as a C _1 function). The resulting extension of the homogeneous solutions, 
and ^J m , written as 

*/m(t, r) = *+ m (t, r) 9[r - r p (t)] + 9j m (t, r) 9[r p {t) - r], (1.2) 

where 9[r — r p {t)\ is the Heaviside function, is a type of weak solution to the inhomogeneous master equation. Thus in 
the gravitational case in RWZ gauge the difficulty with local convergence occurs with the master function itself. We 
show that the use of distributions, or generalized functions (40j . makes possible separate analytic calculation of the 
expected jumps in value and slope of ^em- We further demonstrate that the metric perturbation can be accurately 
numerically computed, including the time dependent magnitudes of delta function terms that appear in some of the 
metric amplitudes in Regge- Wheeler gauge. 

This paper is organized as follows. In Sec. |H] we briefly outline the general mathematical problem of using FD 
techniques to solve for perturbations in the RWZ formalism. We also review the standard parameterization of eccentric 
orbits. Sec. Mil concerns the method of extended homogeneous solutions. We first review BOS's solution for the scalar 
field case. We show then our treatment of more general source terms and extension of the method to gravitational 
perturbations. Sec. IIVI provides numerical results on the computed Fourier-harmonic mode functions, including 
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convergence tests and calculation of radiated gravitational wave energy and angular momentum. In particular, the 
energy and angular momentum fluxes are shown to agree with past published values. More importantly, the method 
is shown to provide a solution to the field and its derivatives that is convergent exponentially fast everywhere. Then 
in Sec. |Vj we show that the equations which allow the metric to be obtained from the master functions, along with 
an understanding of the form of the weak solutions for ^f^ n and ^1^, can be used to determine both the smooth 
and distributional parts of the metric. App. [3] discusses fully evaluated forms of distributional source terms. App. [B] 
gives the details of such source terms for our case of eccentric orbits on Schwarzschild. In App. [C] we concisely 
summarize the metric perturbation formalism in the Regge- Wheeler gauge. We show the construction of gauge- 
invariant master functions of each parity, and provide the spherical harmonic decomposition of the Einstein equations 
and Bianchi identities. App. [D] concludes this paper with a brief discussion of asymptotic expansions used to set 
boundary conditions on the mode functions at large r. 

Throughout this paper we use the sign conventions and notation of Misner, Thorne, and Wheeler [4~i| and use units 
in which c = G = 1. We use Schwarzschild coordinates x^ L — (t, r, 6, ip) except as otherwise indicated. 



II. BACKGROUND ON THE STANDARD RWZ APPROACH TO GRAVITATIONAL 
PERTURBATIONS IN THE FREQUENCY DOMAIN 

In this section we briefly summarize both the standard notation for parameterizing bound orbits on Schwarzschild 
and the usual approach to computing gravitational perturbations using the Regge- Wheeler- Zerilli (RWZ) formalism 
in the frequency domain (FD). The description of the geodesic motion on the background, in terms of various curve 
functions, is used throughout the rest of the paper. The standard FD analysis provides the notation for describing 
the Fourier-harmonic modes, and their normalization, and sets the stage for discussion in Sec. Illll of how gravitational 
perturbations can be returned successfully to the time domain (TD). Here, and throughout this paper, we use a 
subscript p to indicate evaluation along the worldline of the particle. 



A. Bound orbits on a Schwarzschild black hole 



Consider bound timelike geodesic motion around a Schwarzschild black hole (i.e., fi —} 0). We may for the nonce 
use proper time r to parameterize the geodesic, x£(r) = [tp{T),r p (r), p {t), <p p (t)], with the associated four-velocity 
vf* = dx^/dr. On Schwarzschild we take p (r) = tt/2 without loss of generality. The geodesic equations yield 
immediate first integrals and allow the trajectory to be described by the conserved energy £ and angular momentum 
C per unit mass. Alternatively, we can choose the (dimensionless) semi-latus rectum p and the eccentricity e as orbital 
parameters (c.f., [HI, EH)- A third choice would be use of the periapsis and apapsis, r min and r max . We will find all 
of these useful in what follows. The latter two parameter pairs are related to each other by 

^ 2r max r m i n ^ ''max ^min (2 1) 

^(^max ~t~ ^min) ''max ~t~ ''min 

or inversely 

pM pM 



'max — i ' min 

1 — e 1 



(2.2) 



The specific energy and angular momentum are related to p and e by \Tl 

c2 _ (p-2-2e)(p-2 + 2e) r2 _ p 2 M 2 

p(p — 6 — e ) p — 6 — e z 



The geodesic equations provide the following differential equations for the orbital motion and for the time dependence 
of the four-velocity, 



dtp _ t £ dip p _ _ C 



dr f p dr r\ 

where 



^) 2 = K) 2 = £ 2 -c/ p 2 , 



"' = ^ ^?=«* = 3. 3f =(ur = ?-U-, (2.4) 



/(r) = l-— , U 2 (r,L 2 ) = f(l + ^). (2.5) 
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For purposes of numerical integration there is another curve parameter, originally devised by Darwin [42j |. that 
proves useful. Here one introduces a phase angle x that is related to the radial position on the orbit by the Keplerian- 
appearing form 



r pix) 



pM 



1 + e cos x 



(2.6) 



Of course, in the relativistic case x differs from the true anomaly tp. The orbit goes through one radial libration for 
each change Ax = 2tt. The use of x eliminates singularities in the differential equations at the turning points [l5j . 
Note that at x = 0, r p = r m - ln and at x = 77 1 r p — r max- (Also note that in this section we are content with making 
a slight abuse of notation in jumping from r p (r) to r p (x), before ultimately settling on r p (t).) In terms of x the 
equations are 



p 2 M 



dt,. 



d>X (p — 2 — 2ecosx)(l + ecosx) 5 



(P-2) 



4e 2 



p — 6 — 2e cos x 



1/2 



(2.7) 



d<p p _ 

d-X Ip ~ 6 — 2ecosx 



1/2 



and 



Mp 



3/2 



dr„ 



dx (1 + ecosx) 2 



p - 3 - e 2 



-,1/2 



p — 6 — 2e cos x 



We use Eq. (|2.7j) to derive the fundamental frequency and period of radial motion, 



— y > 



(2.8) 



(2.9) 



(2.10) 



It is also of importance to have the average rate at which the azimuthal angle advances, found by averaging the 
angular frequency dip p /dt over a radial libration via 



= JL f f^£p 



Tr 



V dt 



dt. 



(2.11) 



While T r represents the lapse of coordinate time in a radial libration, the time T v — 2tt/Q v has no particular physical 
significance [43j| . Finally, because wave equation source functions contain terms like S[r — r p (t)] and S'[r — r p (t)}, we 
have need of derivatives of r p (t), 



r 2 p (t) - ll 



f 2 

£2 U P> 



r P (t) 



2Mf p f 2 



3M 



C 2 5MC 2 



p p 

where we let a dot signify differentiation with respect to coordinate time. 



(2.12) 



B. The Regge-Wheeler-Zerilli formalism in the frequency domain 



As discussed in the Introduction, we use the RWZ approach to gravitational perturbations and use specifically the 
even-parity Zerilli-Moncrief function ^ c ^ u [13] and the odd-parity Cunningham-Price-Moncrief function [3]. 
See Martel and Poisson [25( for recent discussion and references therein. Both of these functions satisfy wave equations 
of the form 



-W 2 + d^- Mr) 



®e m (t, r) = Se m (t, r), 



(2.13) 



where r* = r + 2M hi(r/2M — 1) is the usual tortoise coordinate. The potential used in Eq. (|2 . 1 3[) is either the Zcrilli 
or Regge- Wheeler potential depending on whether the parity is even or odd, respectively. 

The source terms also depend upon parity but further depend on which specific master functions are chosen. Martel 
and Poisson gave the covariant form of Sf^ n and S£l (see App. [C] for these in Schwarzschild coordinates) that are 
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associated with the Zerilli-Moncrief and Cunningham-Price- Moncrief functions. Martel [3(| derived the detailed form 
of Sp™ n for a point mass in eccentric orbit. Sopuerta and Laguna [35[ derived the detailed form of S^ d for eccentric 
orbits (see also Field et al. [44]). We give in App. |B] detailed expressions for these sources in a form that is useful for 
both mode integrations and metric reconstruction. 

In each case the source term has the following general form 

S im (t,r) = G im (t) S[r - r p (t)] + F tm {t) S'[r - r p (t)}, (2.14) 

where Gi m (t) and F^ m (t) are smooth (differentiable) functions. Note that the source, as written here, differs from 
notation originally used by Martel [3(| (who retained smooth functions of r and t, as in Eq. (11.11) '). Our expression 
uses the delta function, and parts integration, to yield a fully evaluated form along the worldline of the particle (see 
App. [A"| . making Gg m (t) and F im (t) unique functions of time only. 

Eq. (|2.13p can be solved directly in the TD-an approach that has received much attention lately. In this paper we 
are interested instead in extending the reach of FD analysis, and the balance of this section provides a brief review of 
the standard FD solution. We note in passing that a hybrid approach is possible-using FD analysis for low I and m 
modes while using TD calculation for high order modes [45[ . 

On Schwarzschild, eccentric orbits are typically not closed and therefore the motion is not simply periodic as seen 
by an asymptotic static observer. The radial libration is periodic (but not typically sinusoidal) with fundamental 
frequency f2 r . The smooth functions Gi m (t) and Fi m (t), which depend upon the particle's radial and angular mo- 
tion, have terms that are periodic with fundamental frequency Q r , but also involve a term that is proportional to 
exp[— irmp p (t)]. This latter term comes from restricting the spherical harmonics Y^ m (9,(p) with 5[(p — ip p (t)]. The 
function ip p (t) advances with an average rate Q v , but is modulated (in an eccentric orbit) by a function A.ip(t) that 
is periodic with fundamental frequency fl r . Hence, the source Sg m (t,r), and therefore the field ^g m (t,r), can be rep- 
resented by a Fourier series with fundamental frequency f2 r , but multiplied by a phase factor that advances linearly 
with rate il^. These fields would appear simply periodic to an observer whose frame rotates at rate £l v [l5j . To a 
static observer, a given mode i and m will have a spectrum of harmonics offset by mO v ; taken together the full field 
will have a two-fold countably infinite frequency spectrum, 

u) = Lo mn = mVL^ + nVL r , m, n G Z. (2-15) 

Accordingly, the wave equation (|2.13p Fourier transforms into a set of ODE's, 

d 2 



,2 



Rlmn{r) = Z £mn (r), (2-16) 



where Rt mn (r) and Zg mn {r) are Fourier harmonic amplitudes 

R lmn {r) = i- / * dt * /m (t, r) e*"- 4 , Z tmn (r) = i- / ^ ' dt S em (t, r) e*"™*. (2.17) 

Jo J-r JO 

The series representations of ^e m (t,r) and Sg m (t,r) are 

oo oo 

*/ m (t,r)= fltaWe-**"', S tm (t,r)= ^ Z lmn {r) e -™™\ (2.18) 

n=— oo n= — oo 

and are subject to the usual provisos of Fourier theory regarding for what r Eqs. (|2.18[) converge to the original 
functions. 

In order to find the solution to Eq. (|2 . 16[) , we start by solving the homogeneous version of that equation, obtaining 
two independent solutions. Using the terminology of Galt'sov [46| (see also for a clear presentation of basis 
modes), the Rg mn {r) solution is computed by setting a unit normalized "in" wave boundary condition of 

RJ mn (n^-oc)=e-^% (2.19) 

near the horizon. Similarly, the Ki mn {r) solution arises from setting a unit normalized "up" boundary condition of 

RLn(r*^+™) = e^- r ', (2.20) 

at large r*. Formally, these homogeneous solutions are both valid in the entire range 2 M < r < oo. The standard 
method of integrating the Green function and source (the method of variation of parameters) gives the solution to 
the inhomogeneous equation (|2.16|) . 

Rtmn{r) = 4mn( r ) R tmn( r ) + C 7mn ( r ) R l>mn M > ( 2 - 21 ) 
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where 



and 



(r) - 1 r dr , X7 mn (r')Z imn (r>) 1 [ r ™ , RL n (r')Z imn (r') 

l[ ) ~W emn L dr fir') ' Ci ™ [r} ~ W imn I /(r') ' {2 - 22 > 



tit — 6- dR lmn _ + dR imn m oq \ 

M^ m „ = n emn —j- tt emn —j^—, 

is the Wronskian. Outside the source libration region, Eq. (|2.21|) reduces to the normalized homogeneous solutions 
that are properly connected through the source region, 

where C^ nn are the values of cf mn (r) evaluated at the ends of the range of the source, 

Ctmn = C tmn ( r max) , C \ mn = C £mn (r m i n ) ■ (2.25) 



III. THE METHOD OF EXTENDED HOMOGENEOUS SOLUTIONS IN THE GRAVITATIONAL CASE 



A. Brief review of Barack, Ori, and Sago's method of extended homogeneous solutions 



As a model problem, Barack, Ori, and Sago (BOS) considered the scalar field $ produced by a scalar point charge 
in an eccentric orbit on a Schwarzschild background. The spherical harmonic amplitudes (j>i m {t, r) — r$>g m (t, r) of the 
scalar field satisfy RWZ-like equations fully analogous to Eq. (|2.13p but with source functions that only depend upon 
a Dirac delta function, 

S^ = C im (t,r)S[r-r p (t)]- (3-1) 

Here Cf m (t, r) is a smooth function that is derived from the particle's point-like charge density p. 

With a delta function source the amplitudes <fie m (t,r) are left piecewise continuous [C ) at the instantaneous 
particle location r p {t) but lose all differentiability there. BOS argued that this behavior, while surmountable in TD 
calculations, would cause difficulties for Fourier synthesis in FD calculations. As they convincingly demonstrated 
with their first two figures, while 4>£ m (t, r) converges exponentially fast outside the radial libration region, the Gibbs 
phenomenon is responsible for a very slow convergence of ^>e m (t,r) between r m ; n and r max . Furthermore, the radial 
derivative dr<pim is discontinuous at r p (t) and suffers the full effects of the Gibbs phenomenon-the Fourier series 
converges to the mean value at the discontinuity and partial sums (~N < n < N) overshoot in the limit as both 
N — > oo and r — > r p (t) ± . This behavior is a serious obstacle to straightforward use of FD calculations in SF 
regularization. 

As a solution to this problem, BOS developed the method of extended homogeneous solutions (EHS). Their method 
involves using the Fourier-harmonic modes of the homogeneous equation in the FD to synthesize homogeneous solu- 
tions 4'i m (t> r ) ano - < / > foi(^' r ) t° the TD wave equation. The Fourier convergence of these homogeneous solutions is 
exponentially rapid. While these solutions exist in the entire radial domain (2M < r < oo), ordinarily 4>J m (t,r) and 
4>~^ m {t, r) would be viewed as meaningful in their respective source-free regions, r < r m ; n and r > r max . The heart of 
the BOS method lies in extending both of these solutions into the region of radial libration up to the instantaneous 
position of the particle. 

BOS demonstrated the method numerically using the monopole term of $. A key condition for success of the 
method is that, as N — > oo in the partial sums, one finds 

to* <t>Im&r)= lim # m (t,r), (3.2) 

as expected analytically. This was observed numerically and the method as a whole converges rapidly since the FD 
solution of the inhomogeneous equation is never summed. BOS went on to argue that the method could be extended 
to any I and m for scalar, electromagnetic, or gravitational fields. 
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B. Application to gravitational perturbations 



In this section we detail our application of the method to the gravitational case in RWZ gauge. It is worth first 
observing the magnitude of the problem to be circumvented. Given the gravitational source (|2.14[) . and the solution 
to Eq. (|2.16l) afforded by Eq. (|2.21l) , the standard approach would represent the inhomogeneous solution to the master 
equation (|2.13l) by 



tt /m (t,r)~*i£(t,r) 



+N 

£ 

v=-N 



Rimn(r) e 



N^oo, (3.3) 



to indicate that the equality between the actual solution ^ 



and holds almost everywhere 



where we use the 
for N — >■ co. 

Looking ahead somewhat, we use our numerical code to obtain a particular spherical harmonic amplitude, ^22 {t, r) 
(I = 2, m — 2), and its radial derivative, d r &22(t, r). We can also use the code to assemble the standard partial 
Fourier sums (see FIGs. [T]and[2|). We find that the Gibbs problem with the standard approach is significantly worse in 
the gravitational case ( in Regge- Wheeler gauge ) than it is for the scalar field. In the present case the field itself has a 
discontinuity and the radial derivative is both discontinuous as r — ¥ r p (t) and also has a delta function singularity at 
r p (t). The left panels of FIGs.Q]and[3]are familiar; the partial sums have difficulty representing the jump discontinuity 
and overshoot the exact solution (solid curve). In the right panels, the singularity at r p (t) wreaks havoc on the ability 
of the Fourier synthesis to represent the exact solution. 
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FIG. 1: The standard FD approach to reconstructing the TD master function and its r derivative. The left panel shows $22 
and the right shows 8^22 at t = 51.78M for a particle orbiting with p — 7.50478 and e = 0.188917. This figure is analogous 
to FIG. 1 of BOS [3^ | . Partial sums are computed with Eq. (|3.3|l and shown for different N. For contrast we plot the converged 
solution from the new method with a solid curve (see FIG. [3}. The arrow in the right panel gives a notional representation of 
the delta function singularity present in d r ^22\ the amplitude of this singular term is related to the jump in ^22 seen in the 
left panel. 



On a bright note, outside the range of the source, the standard solution converges exponentially fast. Nevertheless, 
in the source region between r m j n and r max the convergence will be algebraic in general and disastrous at the location 
of the particle. A discontinuous (or worse, singular) function cannot be accurately represented by a sum of smooth 
functions. 

We now generalize the EHS method to the gravitational case. We start by recognizing that Rf mn from Eq. (I2.24[) 
are valid solutions to the homogeneous version of Eq. (j2 . 16[) throughout the entire domain outside the black hole, 

4nW = C4n«. r > 2M. (3.4) 
Next, we use these to define the time-domain extended homogeneous solutions, 

= E^nW^""', r > 2M, (3.5) 



9 



2.0 - 



1.5 



1.0 



0.5 



0.0 



1 1 


1 


11 


i i i i i _ 

TV = 10 


- 

- 






TV — 25 ~ 

TV — 50 

*™ s 






// / 

V 

• 




- 


,1 




- 




d 






^ ' . . . • 

i 


1 


iii 


Ill 



6.5 



7.0 



7.5 



8.0 



CD 



10 



12 



II 

\ 
\ 
\ 

1 

- a 


i 

N - 10 

\ TV — 25 - 

! TV — 50 

\ 

I : 


i i i 1 i V i i 1 i i i 1 i 



6.5 



7.0 



7.5 



8.0 



r/M 



FIG. 2: An alternate view of the behavior presented in FIG. [T] A change in the scale in the left panel emphasizes the Gibbs 
overshoots in ^22- On the right, a zoom-out of the vertical scale more clearly indicates the attempt of the Fourier synthesis to 
capture the delta function at r p (t). 



which result from inserting Rf mn into Eq. (|2.18p . The central claim is then that for any t and r the actual solution 
to the inhomogeneous wave equation (|2.13p is given by 

* tm (t, r) = *f HS(i, r ) = tt+Ji, r)9[r- r p (t)] + 9J m (t, r) 9 [r p (t) - r] . (3.6) 

The argument made by BOS can be extended to the gravitational case and goes as follows: 

• We denote the desired true solution of the inhomogeneous wave equation as ^i m . Outside the domain of the 
source (r < r min , r max < r) = = *f„ S because there R tmn = Rf mn - 

• It is assumed that \fV m (t, r) is analytic in the entirety of the two regions 2M < r < r p (t) and r p (t) < r (excluding 
only a neighborhood of r p (t)). 

• Because the homogeneous solutions are expected to be analytic everywhere, ^f^ s (t,r) will be analytic in 
the two regions discussed above (excluding only a neighborhood of r p (t)). (See the extended discussion BOS 
have about this.) 

• Because 'f i m and I J / ™ S are identical outside the region of libration, and they are both analytic everywhere up 
to the location of the source, they must be equal over that entire domain. 

Here we provide an additional justification for the assumed form of the solution given in Eq. p.6[) . The source 
term of the wave equation is a distribution, or generalized function [ifij . Accordingly, any solution of Eq. (|2.13[) 
will be a weak solution-a generalized function itself-with loss of (classic) differentiability at the singular point r p (t). 
To determine the suitability of Eq. (|3.6I) as a solution of Eq. (|2.13l) , we generalize the concept of differentiation to 
encompass distributions. Thus, for example, d9(z)/dz — 5(z). We can then take Eq. p. 61) as an ansatz, substitute in 
Eq. (|2.13p . and determine what conditions are required that it be a (weak) solution. For clarity, in the rest of this 
section we suppress the t and m indices. 

Rather than use the RWZ equation as it stands, we introduce a coordinate transformation to fix the position of the 
singularity. Defining z = r — r p (t), t = t, the derivatives transform as d Tt = f{r)d z and dt = c% — r p d z , and the wave 
equation (|2.13l) becomes 

L(tf) = -flf* + (f 2 - fl) dl^> + 2f p d- t d z ^ + (r p + - V* = G6(z) +F6'(z). (3.7) 

Now we assume that ^> has the form given in Eq. (13.6[) and substitute it into Eq. (13.71) . The functions \E' + and 
are differentiable and satisfy the homogeneous equation, L(5' ± ) = 0. A term of the form L(^ + ) 9{z) + L(ty~) 6{—z) 
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FIG. 3: The EHS approach to reconstructing the TD master function and its radial derivative. As in FIG. [TJ we give ^Z™ 8 
and a r *™ S at t = 51.78M for a particle orbiting with p = 7.50478 and e = 0.188917. Partial sums of *f 2 HS are computed from 
Eq. (|3.5p . with a range of — N < n < N. The full alm its r derivative result from N = 10, which gives agreement in the 

jumps in Vtfl 13 an d d r ^22 to a relative error of 10~ 10 . On the right, the presence of a delta function singularity is notionally 
depicted with an arrow. The time dependent amplitude of this singularity is separately computable from the jump in ^22- 



appears in (|3.7[) and drops out. Other singular terms remain, created by derivatives of the Hcavisidc function, and 
we are left with 

(/ 2 - r 2 p ) (ld r *} p S(z) + m p 6\zj) + 2f p $([*] p *(*)) + (r p + (fd z f) ) [#]„ S(z) = GS(z) + F S'(z). (3.8) 
where 

M P (t) = (t,r p (t)) - *" (t,r p (t)) , ld r *l P (t) = d r V+ (t,r p (t)) - d r *~ (t,r p (t)) (3.9) 

are the jumps in VP and d r ^ at z = 0. Naively, we might expect that we can simply equate the coefficients of 6 on 
the two sides of Eq. 13.81 while doing the same with the 6' coefficients. However, the 6' term on the left hand side 
must first be fully evaluated (as a function of time) at the location of the particle. To do this, we use the identities 
in Eqs. (|A1|) and (|A5|) . which leaves 

(f P - r 2 p ) M P 8(z) + { f p - +1) m P S'(z) - 2 (f p d z f p ) m p 8(z) + 2f p 9 t -([*] p ) 5(z) 

'r P + (fpd z f P ) ) Mp S(x) = G S(z) + F 5'(z), (3.10) 



(' 



where f p = f{r p (t)). Note that there is no comparable expansion on the right side from the FS'(z) term because 
F is already fully evaluated at r = r p (t), by design. From here, we read off the jumps in ^> and its r derivative at 
r p (t) from the coefficients of §' and 5, respectively. Returning to Schwarzschild coordinates and using Eqs. (|2.12l) to 
remove r„ and rz terms, we find 



v 



F , ,.x 



J p ^ p J p ^ p 



(3.11) 



From the standpoint of the original coordinates, the partial time derivative <9j becomes the convective, or total, time 
derivative along the particle worldline. 

These jump conditions amount to internal boundary conditions that are necessary conditions on a solution to the 
inhomogeneous wave equation in the TD. They were discussed by Sopuerta and Laguna [35| and also later, with 
corrections, by Field et al. [44|. In our FD-based calculations, they provide a powerful check on our transformation 
of the solutions back to the TD. Given the indirect way in which the Fourier transform of the source Si m determines 
the Fourier coefficients of the extended homogeneous solutions, considerable credence is lent to the method in seeing 
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the partial sums of converge toward satisfying these jump conditions. Secondarily, the jump conditions provide 

useful stopping criteria in the numerical method (see Sec. II V Cj) . 

While not a focus of this paper, we consider briefly TD simulations. There, to find a unique solution the internal 
boundary conditions must be augmented with initial data on a Cauchy surface and, potentially, outer boundary 
conditions. Care must be exercised to switch on the source smoothly in the (near) future of the initial value surface 
plil ] (also Lau, private communication). Additionally, imposed initial data will not typically match long term periodic 
behavior induced by the source, and transients will sweep through the system for several dynamical times. In contrast, 
in the FD approach, the proper outgoing and downgoing behavior at the outer boundaries is built in from the outset 
and only the steady state, periodic behavior is obtained. 



C. Computing normalization coefficients in the gravitational case 

Finally, we provide some details on how the singular source is integrated to provide the matching normalization 
coefficients Cf mn and Cj mn that are used in Eq. p.4[) . BOS detail the calculation of normalization coefficients for the 
scalar monopole in their App. C. The gravitational case follows the same general idea, but involves some technical 
differences and challenges. We start by combining Eqs. (|2.25j) and (|2.22|) . giving 

^^r ^y) , (3 . 12) 

The FD source term Zi mn (r) comes from plugging Eq. (|2.14[) into Eq. (|2.17[) . yielding 

Z lmn (r) = i- dt (G lm {t) 5[r - r p (t)} + F (m {t) S'[r - r p (t)])e iw ™"*. (3.13) 

The equivalent integral BOS present for the scalar monopole is their Eq. (C2), which they evaluate immediately by 
changing the integration variable from t to r p . Here, with a derivative-of-the-delta function present (in RWZ gauge), 
the immediate evaluation of this integral produces terms that are singular at the turning points (f p — 0). These terms 
are no problem analytically, but they are troublesome when performing the final numerical integration. We therefore 
find it is advantageous to delay this integration. Plugging our expression for Zi mn in above, we have 

Cf m n = ~uf~~~~T~ f ^ dr %^ F dt (G em (t) S[r - r p (t)] + F tm (t) S'[r - r p (t)\) e**--'. (3.14) 

Wemn-l-r J r min J\ r ) JO V ' 

In order to avoid the singularity at the turning points, we switch the order of integration. The integration of the 
delta function itself is then straightforward. The derivative of S term requires an integration by parts. Because of the 
compact support of the source term, we can extend the range of integration and no surface terms appear. We are left 
with 



a 



± 



1 



ran tj/ rp 

VVimn-i-r JO 



jRj mn irv)G im {t) + (^- 2 Rj mn (r P ) - J p dRfm d f p) ) Fimit) 



dt, (3.15) 



where we use a p subscript to indicate evaluation of a quantity at r = r p (t). Our final integral is analogous to Eq. (C7) 
in BOS. 

Here is a summary of key details of the application of the method in the gravitational case: 

• The EHS method, applied to the gravitational case, gives exponentially converging solutions to Eq. (|2.13[) 
everywhere, including the location of the particle. (See FIG. HI) 

• Working in Regge- Wheeler gauge, the gravitational TD source term contains a delta function and a derivative- 
of-the-delta function, which cause ^e m to exhibit a jump and d r ^i m to exhibit both a jump and a delta function 
singularity at the particle's location. (See FIG. 03) In the scalar case, the field is piecewise continuous at the 
particle, with a jump in the r derivative. (See FIG. 3 in BOS.) 

• Eq. (|3.15[) is valid for all radiative multipoles (£ > 2). The t = 0, 1, modes must be handled separately. 

• Martel's [3Cj Ge m (t,r) and Fg m (t,r) from Eq. (jl.ll) are not in fully evaluated form. As discussed in App. [5] for 
a given multipole, unique functions of time F^ m (t) = F lm (t, r p (t)) and Gt m (t) = Gi m (t, r p (t j) — d r F e . m (t, r p {t)) 
emerge after fully applying the delta function constraint. We use the tilde to distinguish fully evaluated coeffi- 
cients. 
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• In practice, we take advantage of the fact that some of the functions in the integrand of Eq. (|3.15l) are even 
over the period of radial libration, while others are odd. Then, rather than integrating over t from — > T r , we 
can limit the range of integration to — > T r /2. Further, we change variables to x^ as shown in Sec. Ill Al and 
integrate from — > n. 

• For V&fl? 11 we use the Zerilli-Moncrief master function, and for \E'°^ l d we use the Cunningham-Price-Moncrief 
master function. This formulation works for any master function that obeys a Regge- Wheeler-like equation and 
has a source term that can be written in the form of Eq. (|2.14p . 



IV. NUMERICAL METHOD AND RESULTS FROM MODE INTEGRATIONS 



A. Algorithmic roadmap 

Here, we explain the specific steps our code takes to solve the inhomogeneous wave equation (|2.13[) . There are 
several stages to the process, and at each step we compute at least one more order of magnitude accuracy than is 
needed at the subsequent step. The code is written in C, and we use the Numerical Recipes adaptive step size fourth 
order Runge-Kutta integrator (48j . 

1. Specify an orbit through a choice of the semi-latus rectum p and eccentricity e. 

2. Numerically integrate Eqs. (|2.10l) and (|2.11l) to get the fundamental frequencies of the system, f2 r and ft v , and 
hence uj mn = mVL v + n£l r . 

3. Choose a specific £ and m. If I + m is even (odd), use even (odd) parity potential and source terms. Choose 
starting n. (See Sec. HVCl ) 

4. Solve the homogeneous version of Eq. (|2.16[) to get unit normalized radial mode functions, R^ mn , in the source- 
free region: 

• Use the asymptotic expansion (see App.[D]) to set an "up" plane wave boundary condition at ?'» — > +oo, as 

in Eq. (|2.20[) . Numerically integrate up to the region of the source at r" lax to get R\ mn - (We let r " lm / max 
be the r* value corresponding to r min / max .) 

• Use a convergent Taylor expansion to set an "in" plane wave boundary condition (Eq. (|2.19l) ) at modestly 
negative r*. Numerically integrate up to the region of the source at r™ m to get Rj mn - 

5. Solve the homogeneous version of Eq. (|2.16j) to continue the unit normalized radial mode functions, Rf mn , into 
the source region, while also computing the normalization coefficients C^ mn : 

• Simultaneously integrate Eqs. (|2.16[) and (|3.15[) from \ = — > it (cquivalently t = — > T r /2 and r = 
fmin -> r malt ). This gives RJ mn in the region of the source and C/ mn . 

• Simultaneously integrate Eqs. (|2.16[) and (|3.15[) from \ = —n — > (equivalently t = —T r /2 — > and 
r = r max r min ). This gives Rj mn in the region of the source and Cj mn . 

As discussed in Sec. IIII C[ the integrand in Eq. (|3.15l) contains parts which are even and parts which are odd 
over the radial period. By keeping the correct terms, we can get away with efficiently integrating over only half 
the period. 

6. Use the coefficients to normalize the homogeneous solutions outside and inside the range of the source, as in 
Eq. (EU). 

7. Assess whether there is convergence of the partial sum over n. (Again, see Sec. IIV CI ) 

• If yes, we are finished with this £, m mode. 

• If no, return to Step U with the next n. 
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B. Energy and angular momentum fluxes at r* = ±00 



To evaluate the energy and angular momentum fluxes at r* = ±00 we use the Isaacson stress-energy tensor. The 
energy and angular momentum fluxes, for each £, m mode, can be written as 49] 



1 (1 + 2)! 
64tt (£-2)\ 



im (£ + 2)1 
64tt (£-2)\ 



*f m (*»0*L*(<,r). 



(4.1) 



Here, an asterisk signifies complex conjugation. (We use VE^™ 11 when £ + m is even and when £ + m is odd. In 
general there would be contributions from both x I'^™ n and for each mode, but our choice of 9 P = n/2 leads to 

one of these functions vanishing for each £ and m combination.) In terms of FD amplitudes the expressions become 



eL 



£< 



1 (1 + 2)! 
64tt (£-2)! 

m (1 + 2)! x ^ 

54^ (^ — 2)! L "" m ~ Il * nn '"* Tm 



1 0^.^.1 1> 



(4.2) 



As is well known, the fluxes must be suitably averaged over time or space to obtain meaningful, invariant results. We 
average these quantities in time over one radial oscillation, which yields 



E 



1 



2)! 



64tt (£-2)! 



r± r± 



f ± \_ m (£ + 2)1 



LO r , 



r± r± 



(4.3) 



Here, we have also introduced R lmn = C hnn R lmn 



Re 



fluxes at r* = ±00, leaving 

,±oc\_ 1 (* + 2)! 



where J Ln( r 



± 
Ci 

1 as r 



/ 64tt(£-2)!^' 

As discussed in App. [Dl we can write the radial function as 
±00. Therefore, if we set J^„„ = 1, we can evaluate the 



El 



64tt (£-2)1 



E 



mra | ^ Imn \ 



±00 



m (^ + 2)! 



64tt (£ - 2) 



c, 



± 



(4.4) 



C. Code validation 



To compute the total energy and angular momentum fluxes, we must sum Eqs. (14.41) over £ and m. The resulting 
expressions are formally over the ranges 2 < £ < 00, — £ < m < £, —00 < n < 00. When computing _E and 1/ 
numerically, we put limits on each of these sums. To begin with, the low £ modes matter more than the high ones. 
But, the more eccentric an orbit, the more £'s must be computed to achieve the same precision in our final values. 
For the orbits we considered in Table HI in order to achieve a relative precision of 10 -12 in our final flux values, the 
highest £ necessary was £ — 29. (See Sec. IIVD1 ) 

Because of the symmetry of the spherical harmonics, the fluxes from any given — m mode are equal to those from the 
corresponding +m mode. Therefore, we fold the negative m modes over onto the positive ones, and simply multiply 
each positive m mode by two. Additionally, as £ gets larger, it is no longer necessary to compute all m values. As can 
be seen in Table IIII1 for a given £, the largest E^ H and LuI H contributions come from the m = £ mode. We start 
at m — £ and decrement m until the fluxes are no longer significant. For low £ values we still wind up computing all 
< m < £, but as £ increases, we need progressively fewer m modes. 

Determining the necessary n's is a bit more involved. For a given £ and m, there is a range, n min to n max , over 
which we sum in order to achieve our desired precision. Looking at Table Hill it is evident that when m = 0, the range 
of n is essentially centered on 0. For these modes, we start with n = 0, and compute fluxes for all positive modes. 
When we have seen no change to any of the flux values (at a pre-specified level of precision) for several consecutive 
modes, we stop and repeat the process for the negative n's. As m increases, this range of n's shifts more and more 
into the positive. For any £, the m — £ mode has far more positive n modes than negative. Eventually, £ becomes so 
large that n m j n > for the m = £ mode. For modes where we suspect that n m i n > 0, we find it advantageous to start 
with a rough sweep of a large range of possible n values. We calculate E^ nn (the energy flux at r = +00 from one 
n mode) to low precision for a small number of n, spaced out over this range. The n for which we find the largest 
Eimn will be near the center of the n m ; n to n max range. We then perform our high precision mode integrations for all 
significant n values above and below this n. 
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If we are interested in a local calculation (as one would perform for a SF evaluation), we have a different method for 
determining which n's are significant. We still use the energy fluxes to find the approximate center of the significant 
n range, but for the "breaking condition" we compute n's until the jumps in and d r ^i m converge properly, as 
follows: 

• Use Eq. (|3.5p to compute a partial mode sum approximation of both ^f m (t,r p ) and d r ^ff m (t, r p ) for a large 
number of times tk throughout the orbit. 

• Numerically evaluate the jumps in those partial sums 

[*Ll p = *tm (*> r p ) - 9J m (t, r p ) , fdr^J p EE d r *j m (t, r p ) - dr*J m (t, T p ) , (4.5) 
for those times tk- 

• Compute the analytical values of and [<9 r v E'^]] p derived in Sec. MI 51 for those times tk- 

• If [*<m] = I^fonlp ano - I^^imlj, = lA^^m] at au times tk, to a chosen precision, we have computed 
enough n modes. 

• Otherwise more n modes are needed. As in the flux computation case above, we perform the mode calculations 
for the n values above our starting n, and once that partial sum has converged to our desired precision, we solve 
for the n's below our starting n until the jump values agree. 



D. Results 



One of our most important results is the exponential convergence of ^f^ s and its r derivative at the location of 
the particle. FIG. [3] shows a partial sum of these two quantities converging after only a few modes. Compare this to 
FIGs. Q]and[2j which shows the standard FD approach. In particular, note in those figures the failure of the standard 
approach to accurately represent d r ^g mi even after a large number of modes. This function is particularly badly 
behaved in the standard approach as smooth functions attempt to capture a delta function. 

Also of note is FIG. [4] which shows that the convergence from the method of extended homogeneous solutions is 
indeed exponential, all the way up to the location of the particle. Fast and accurate computation of ^ i m and d r ^£ m 
at r p (t) will eventually be critical for reliable local SF calculations. 
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FIG. 4: A plot of the convergence of the master function using the two methods. For a particle orbiting with p = 7.50478 and 
e = 0.188917 at t = 51.78M we compute the master function $22(n m ax) by summing over modes ranging from — n max < n < 
n m ax for n max = 15. We plot the log of the difference between ^22(n max ) and the partial sum ^22(iV), for different N < n max . 
For the standard approach (left), we see exponential convergence in the homogeneous region, but only algebraic convergence 
in the region of the source. The method of extended homogeneous solutions (right) yields exponentially converging results 
at all points outside and inside the region of the source. The method of extended homogeneous solutions gives exponential 
convergence for <9 r 'I y ™ S as well. 
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4~ (m 2 /m 2 ) 


Ef m (M 2 /n 2 ) 


LTm (M/fl 2 ) 


L? m (M/fj, 2 ) 


p = 7.50478, e = 0.188917 
This Paper, £ max = 23 
Fujita et al. 


3.16899989185 xlO^ 4 
3.16899989184 xlO -4 


5.23247295625 xlO" 7 
N/A 


5.96755215609 xlO" 3 
5.96755215608 xlO -3 


8.71943028067 xlO' 6 
N/A 


p = 8.75455, e = 0.764124 
This Paper, i max = 29 
Fujita et al. 


2.12360313360 xlO" 4 
2.12360313326 xlO" 4 


2.27177440621 xlO -6 
N/A 


2.77735939025 xlO" 3 
2.77735938996 xlO -3 


2.22781961809 xlO" 5 
N/A 



TABLE I: Total energy and angular momentum fluxes for eccentric orbits, compared with those from Fujita et al., published 








eL (m 2 /m 2 ) 


Lf m (M/ M 2 ) 


Lf m (M/ M 2 ) 


p= 10, e = 0.1 
This Paper 
Fujita et al. 


6.31752474718 xl0~ 5 
6.31752474720 xlO" 5 


1.53365819446 xlO" 8 
1.53365819445 xlO" 8 


1.95274165241 xlO -3 
N/A 


4.48832141611 xlO -7 
N/A 


p = 10, e = 0.5 
This Paper 
Fujita et al. 


9.27335011599 xlO" 5 
9.27335011503 xlO -5 


1.41298859263 xlO~ 7 
1.41298859260 xlO" 7 


1.97465149446 xl0~ 3 
N/A 


2.15617302381 xl0~ 6 
N/A 


p = 10, e = 0.7 
This Paper 
Fujita et al. 


9.46979134556 xlO" 5 
9.46979134409 xlO" 5 


3.55415030147 xlO" 7 
3.55415030114 xlO" 7 


1.63064691133 xl0~ 3 
N/A 


4.20771917244 xl0~ 6 
N/A 


p = 10, e = 0.9 
This Paper 
Fujita et al. 


4.194264692 xlO -5 
4.19426469206 xlO -5 


3.652142848 xlO" 7 
3.65214284306 xlO -7 


5.982866119 xlO" 3 
N/A 


3.518978461 xlO" 6 
N/A 



TABLE II: Energy and angular momentum fluxes for eccentric orbits, compared with those from Fujita et al. [EJ. Partial 
sums for all four orbits are truncated at ^ max = 20 for both papers. Fujita et al. obtained their numbers from integrating the 
Teukolsky equation. We include this table to show the agreement of our horizon energy flux values. 

In order to check our code's accuracy, we computed energy and angular momentum fluxes for circular and eccentric 
orbits. Our circular orbit fluxes agree, mode- by- mode, with published results (e.g. Cutler et al. [50]) to high precision. 
For eccentric orbits, we are only aware that total energy and angular momentum fluxes have been published. Our 
FD results agree with the fluxes at r — > oo of Fujita et al., published in [l4| to at least 10~ 9 . These are included in 
Table HI Fujita et al. have also published horizon energy fluxes [5l|, which we agree with, to at least 10 -9 for a range 
of eccentricities. These are given in Table |TTJ 

For those wishing to reproduce our results, in Table Hill we give mode-by-mode fluxes up to t = 5 at r = oo and 
down the black hole at r = 2M for a particle in orbit with p = 8.75455 and e = 0.764124. Included are the ranges of 
n modes summed over to achieve these results. 

As expected, our code is more efficient for low eccentricities. The first orbit in Table|T](p = 7.50478, e = 0.188917), 
runs in under a half hour on a single processor machine, giving the total flux for all 2 < £ < 23 (although note the 
limits on m and n mentioned in the previous subsection) to a fractional error of 10~ 12 . As e increase, though, run 
times increase greatly. The second orbit in that table (p — 8.75455, e = 0.764124) takes six hours to achieve the same 
accuracy for all necessary 2 < £ < 29. And, when e = 0.9 for 2 < £ < 20 in the last row of Table [TTl we had to raise 
our fractional error to 10~ 10 in order to get a run time of eighteen hours. 

Clearly, as e gets close to 1, FD methods will lose out to TD codes, which handle high eccentricities with more 
ease. Still for < e < 0.9, our run times are not unreasonable when considering the high accuracy we achieve. 

V. RECONSTRUCTION OF THE METRIC PERTURBATION AMPLITUDES 

The full benefit of having complete and highly converged solutions for the master functions lies in using them 
to reconstruct the metric. Ultimately, one wants to use the information, along with an appropriate regularization 
scheme, to compute the self force. A developed approach to doing this is the mode-sum regularization method [EH, 
which makes use of Lorenz gauge. Here we use the information encoded in the master functions to compute accurately 
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t m 




Ef m (M 2 /^ 2 ) 


LTm (M/ fJ 2 ) 


Lf m (M/ M 2 ) 


^min 




2 





1.27486196317 


xlO" 


8 


1.66171571270 


xlO" 


8 
















-74 


76 




1 


1.15338054092 


xlO" 


6 


3.08063328605 


xlO" 


7 


1.44066000650 


xlO" 


5 


2.77518962557 


xlO" 


6 


-62 


78 




2 


1.55967717209 


xlO" 


4 


1.84497995136 


xlO" 


6 


2.07778922470 


xlO" 


3 


1.85014840343 


xlO" 


5 


-47 


82 


3 





2.53527063853 


xlO" 


11 


1.23159713946 


xlO" 


10 
















-84 


85 




1 


9.66848921204 


xlO" 


10 


2.47099909183 


xlO" 


9 


1.93528074730 


xlO" 


8 


2.10622579957 


xlO" 


8 


-66 


87 




2 


6.17859627641 


xlO" 


7 


1.29412677182 


xlO" 


8 


7.54192378736 


xlO" 


6 


1.23105502765 


xlO" 


7 


-48 


93 




3 


3.71507683858 


xlO" 


5 


8.07017762262 


xlO" 


8 


4.67102471030 


xlO" 


4 


7.99808068724 


xlO" 


7 


-34 


99 


4 





1.14820411420 


xlO" 


12 


1.50591364139 


xlO" 


12 
















-80 


80 




1 


4.58377338924 


xlO" 


12 


2.04365875527 


xlO" 


11 


4.50183584238 


xlO" 


11 


1.63314060565 


xlO" 


10 


-77 


94 




2 


1.59253324588 


xlO" 


9 


1.62313574547 


xlO" 


10 


2.40079049220 


xlO" 


8 


1.51029853345 


xlO" 


9 


-51 


93 




3 


2.44084848389 


xlO" 


7 


6.50157912447 


xlO" 


10 


2.91633622588 


xlO" 


6 


6.23354901544 


xlO" 


9 


-34 


106 




4 


1.12530626433 


xlO" 


5 


4.66621235553 


xlO" 


9 


1.37037638198 


xlO" 


4 


4.59633939401 


xlO" 


8 


-31 


114 


5 





2.93546198223 


xlO" 


15 


1.68762144246 


xlO" 


14 
















-94 


94 




1 


1.66341467681 


xlO" 


13 


2.73842758121 


xlO" 


13 


1.99357707469 


xlO" 


12 


2.10992243393 


xlO" 


12 


-77 


92 




2 


1.72172010497 


xlO" 


12 


1.49178217605 


xlO" 


12 


2.59625132235 


xlO" 


11 


1.34307539131 


xlO" 


11 


-63 


100 




3 


1.73935003471 


xlO" 


9 


1.01021973779 


xlO" 


11 


2.26258058740 


xlO" 


8 


9.56603438989 


xlO" 


11 


-46 


109 




4 


9.01787571564 


xlO" 


8 


3.64807139949 


xlO" 


11 


1.06079902733 


xlO" 


6 


3.50623060712 


xlO" 


10 


-29 


121 




5 


3.74854353561 


xlO" 


6 


3.02291684853 


xlO" 


10 


4.47051998131 


xlO" 


5 


2.96568439531 


xlO" 


9 


-19 


130 


Total 


2.10242675876 


xlO" 


4 


2.27174892328 


xlO" 


6 


2.75262625234 


xlO" 


3 


2.22779475534 


xlO" 


5 





TABLE III: Energy and angular momentum fluxes for an eccentric orbit with p = 8.75455, e = 0.764124. Note that we have 
folded the negative m modes onto the corresponding positive m modes and doubled the flux values in this table for m > 0. 



the spherical harmonic amplitudes of the metric perturbation in Regge- Wheeler gauge. The ability to determine the 
metric at all locations, including at the particle location, should serve as a useful starting point for computing the 
SF, either via a gauge transformation or an alternative regularization technique. 

We summarize the metric perturbation (MP) formalism in App. [Cj where the definitions of the master functions, 
$™ and 5 , £^ l d , are given in terms of spherical harmonic amplitudes of the metric and their radial derivatives. We re- 
serve for this section giving the equations, (|5.5[) and (|5.15[) . for reconstructing the metric amplitudes in Regge- Wheeler 
gauge from the master functions. These equations involve first derivatives, and in some cases second derivatives, of 
the master functions. They also involve spherical harmonic projections of the stress-energy tensor. Based on the form 
(|1.2p anticipated in a master function, both of the abovementioned facts contribute to an expectation that the MP 
amplitudes might have point-singular behavior at r p (t) in the form of both S and 8' terms. We show that all potential 
S' terms cancel out. However, in general a MP amplitude might have a functional form 

M(t, r) = M + {t, r) 9{z) + Ar(t, r) 6{-z) + M s {t) 5{z), z = r- r p (t), (5.1) 

where M + (M~) represents a smooth function in the region r > r p (r < r p ), and M is a smooth function of t alone, 
giving the magnitude of the singularity. We examine M s for all six non-zero MP amplitudes in the Regge- Wheeler 
gauge, and find three such terms to be nonvanishing. Throughout the rest of this section we again suppress spherical 
harmonic labels I and m. 

As mentioned the metric reconstruction equations, of each parity, require spherical harmonic projections of the 
stress-energy tensor. For a particle of mass n, traveling on a geodesic of the background spacetime, with four- velocity 
it M , it is 

T»" (x a ) =fi f -^/(t)/(t) S 4 [x - x p (t)} . (5.2) 



In Schwarzschild coordinates the determinant of the metric is g = — r 4 sin 9. After changing the variable of integration 
to coordinate time t, we have 

T " v {xa) = u *(t) IJy s[r ~ rp{t)] 6[lf ~ Mt)] s[e ~ * l2] - (5 - 3) 

Spherical harmonic projections of appear as source terms in the decomposed Einstein equations (App. [C]) and 
these are in turn combined to produce the source terms for the master equations (App. IB]), In the subsections that 
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follow, we evaluate the time dependence of all of the stress-energy tensor projections. We use the definitions 

A(r) = A + — , As e + g)('-l). 
r 2 



A. Even parity 



The even parity MP amplitudes are expressed in terms of $ even and the source terms by (see [301 ]) 

r 2 f 2 



K(t,r) = fd r V even + A% 
"A + l 



h rr (t,r) 



A_ 



-K 



(A + 1)A 



Q 



f 



(5.5) 



fetr(*, = rd t d r ty cvcn + rB d t ^ cvcn 
h tt (t,r) = fkrr + fQK 



A + l 



Q' 



A 



where 



A(A 



1 + A + 



A 1 



3M 
r 



3M 2 



(5.6) 



These equations result from the definition (|C6j) of VP ev0 n and its substitution into the even-parity field equations (IC3|) . 
The even-parity projections of the stress-energy tensor that appear in the equations above are defined by Eqs. (|C4D . 
By enforcing the delta function constraints, they can be written in fully evaluated form (see App.|B|), with each having 
a time dependent magnitude multiplying the radial delta function 



T\t, r) ee q a \t) S[r - r p (t)], Q Q (t, r) = q a (t) 6[r - r p (t)}, 
Q b (i, r) ee q \t) 6[r - r p (t)], Q»(i,r) ee q*(t) S[r - r p (t)], 



(5.7) 



where we use a lowercase q as the base symbol of the corresponding magnitude. With Eq. (|2.4p giving the four- velocity 
it**, the stress-energy tensor and Eqs. ()C4[) can be used to find 



£ 



q u (t) = q rr (t) = 8vr M A (f3 _ y^ 

r plp ^ r p 



£r 



q tr (t) = Snfi-Y* 



167T/J £ 



*(* + l)rf 



2 ¥> 



<f(*) 



167T/X Cf p 

£{£ + l)£rj v ' 



(5.8) 



^(t) = 87T/i -p-^fy*, q s (t)=32nfji 
£ rz 



(£-2)\C 2 f 1 



Py* 



{£ + 2)\ £ r 2 



Here, Y, Y Vl and Y vv are shorthand for the even-parity scalar, vector, and tensor spherical harmonics, respectively, 
evaluated along the worldline at 9 = it/2 and ip — ip p (t). 

Now consider the reconstruction of the MP amplitude K, given in Eq. (15.5[) . Using the expected functional form of 
^ given in Eq. (|1.2I) . K obviously does fit the general form (|5.ip claimed above. In fact, we find 



K^t, r) = fd r *± + A**, K s (t) = f p m p 



r 2 f 2 

pJp 



(A + 1)A 



-q tt = 0, 



(5.9) 



where the vanishing of K s follows from use of Eq. (|3.11[) for \%f\ p , and q tt from Eq. (|5.8I) . Therefore, we see that the 
even- parity metric function K in Regge- Wheeler gauge is (only) a C _1 function at the location of the particle. 
Using the same approach to evaluate h rr in Eq. (|5.5I) we have 



h tr(^ r ) = J2 



A + l 



A,, 



h S rr (t) = -flKlp = r P ld r n P + -"Mp 



fp 



(5.10) 



Here, we have extended in a natural way the use of the [ ] p notation to let fKj p represent the jump in K at z = 0. We 
find that the Regge- Wheeler metric function h rr is not only discontinuous across r p (t) but also has a point-singular 
term, which is an artifact of Regge- Wheeler gauge. 
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The htr function is more subtle than the previous two. Looking at Eq. (I5.5|) . we need the following terms involving 



-B deil = rB d t ^ + 9(z) + rB dt^~ 9{-z) ~ r p B p r p S(z), 
rd t d r V = rd t d r q + 9(z) + rd t d r ^-Q{-z) + 



(5.11) 

5(z)~r p r p M P S'(z). 



On the right side of these equations we have evaluated all the 5 and 6' coefficients at z = with Eqs. (IA1|) and (|A5|) 
(fully evaluated form). The singular terms that arise in these expressions can be grouped with the similarly singular 
contributions from the source terms, 

r 2 r 2 
A + l^ \ + i q {h 

3 , (5.12) 



r3f -8 t Q« - 1 



(A + 1)A ^ (A + 1)A P 



, dq u 3Ar„ + 12Mr p — 4AMr p — 18M J .. 

+ — aZ ™ 



^-(ATik^^- 



Upon carefully checking the time dependence of q u and the jump in ^P, we find that the 5' terms cancel out. There 
are multiple 6 terms, but after using the expressions for l^fj p , [9r*]j9 in (|3.11l) and the relevant g's in (|5.8I) . most of 
the terms cancel and we are left with 

ht r (t,r)=rd t d r V ± +rBd t y ± , hf r (t)=£ 2 -^qK (5.13) 

Finally, the hu term is simple. We insert Eq. (|5.10p into the field equation for h u and get 

hf t (t, r) = fh± , hg(t) = f 2 h s rr + f p qK (5.14) 

So, we see that in Regge- Wheeler gauge K is C _1 with no singularity along the worldline of the particle, but the 
three even-parity MP amplitudes in the "t, r sector" have point-singular artifacts given by Eqs. (|5.10l) . (|5.13p . (|5.14[l . 

B. Odd parity 

Once \l/odd has been computed, the odd-parity MP amplitudes can be reconstructed via 

ht(t, r) = ±d r (r* odd ) - T -J-P\ h r (t, r) = ^jd t ^ odd + ^jP r , (5.15) 

(see [3ljj ) . These equations follow from the definition (|C14I) and its substitution into the odd-parity field equations 
(|C11[) . Similar to before, we define the lowercase p's to be the time-dependent magnitudes of the radial delta function 
after fully evaluating the odd-parity projections of the stress-energy tensor 

P a (t,r)=p a (t)S[r-r p (t)}, P(t,r) = p(t) S[r - r p (t)]. (5.16) 

Also as before, we use the time dependence of the four-velocity and the stress-energy tensor to determine these 
magnitudes for eccentric motion on Schwarzschild, 

t 16nn £ 16nfi C f p {l-2)\C 2 f p 

Here, X v and X vv are shorthand for the odd-parity vector and tensor spherical harmonics evaluated along the 
worldline at 9 = n/2 and if — ip p (t). 

Now, as in the even-parity case we can analyze the local structure of the MP amplitudes. We again assume \1/ to have 
the form Eq. (jl.2[) . Plugging the relevant expressions into Eq. (15 . 15[) for the odd-parity MP amplitude reconstruction, 
we find that all the point-singular parts cancel out exactly, leaving 

hf(t,r) = £d r (r* ± ), h s t (t) = 0, 

2 r ( 5 - 18 ) 
h${t,r) = —d t * ± , h s r (t)=0. 

So, we see that the odd-parity MP functions in Regge- Wheeler gauge are smooth as they approach r p (t) with only a 
finite jump at that point. 

FIG.[S]summarizes these findings graphically, for both even and odd parity, using several specific spherical harmonic 
modes. 
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FIG. 5: The EHS approach to reconstructing the TD MP amplitudes. We consider a particle orbiting with p — 7.50478 and 
e = 0.188917 at t = 80.62M. The left plot shows the odd-parity MP amplitudes /if 1 and /if 1 . The right shows the even-parity 
ft.f|, /iff, /iff, and K 22 . Note that the amplitudes /iff, /iff, and /iff are singular along the particle's worldline, as indicated by 
arrows in the plot on the right. The magnitude of these singularities are given in Eqs. (|5.10|) . (|5.13p . ()5. 14f) . The remaining 
three MP amplitudes approach the particle location smoothly, and have only a finite jump at that point. 



VI. CONCLUSION 



We have achieved two main results with this paper. First, we have shown successful application of the method of 
extended homogeneous solutions to gravitational perturbations from a small mass in eccentric orbit about a massive 
Schwarzschild black hole. In doing so, we accurately computed the master functions in the Regge-Wheeler-Zerilli 
formalism in the frequency domain and transformed these fields back to the time domain. With this method we 
achieved exponential convergence of the master functions and their derivatives for all r including the instantaneous 
particle location r — r p (t). 

Our second important result is the reconstruction of the metric perturbation amplitudes in Regge- Wheeler gauge 
for arbitrary radiative modes. In addition to computing the smooth parts of these amplitudes, we have derived the 
time dependent magnitudes of point-singular terms that reside at r p (t) in some components of the metric. This full 
and accurate knowledge of the spherical harmonic amplitudes of the metric at, and near, r p (t) lays the groundwork 
for one or more subsequent approaches to local computation of the self-force. 
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Appendix A: The fully evaluated form of distributional source terms 



In the RWZ formalism for perturbations generated by an orbiting point mass, the master equations have distri- 
butional sources with both delta function and derivative-of-delta function terms. Reduced by spherical harmonic 
decomposition, these distributions have support only along a one-dimensional timelike worldline r — r p (t) within a 
two dimensional domain. The delta function's behavior is still elementary, 



a(t,r)6[r-r p {t)]=a(t,r p (t)) 6[r - r p (t)] = a(t) S[r - r p {t)}, 



(Al) 
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where a(t, r) is assumed to be a smooth function and we use the notation a(t) to indicate the one-dimensional function 
that results from restricting (or fully evaluating) a(t, r) with the delta function. At any stage in a calculation a delta 
function can be used to fully evaluate all smooth functions that multiply it. Under an integral the result is obvious 



a(t, r) 8[r — r p (t)] dr = a(t), 



(A2) 



with the resulting function of time being unique. Occasionally, there is need to differentiate such a function. The 
total derivative is related to derivatives of the original function by 



da 
~dt 



d t a(t, r) + r p d r a(t, r) 



r=r p {t) 



where on the right hand side we differentiate first and evaluate second. 

Of more interest is the behavior of 8' [40j. Differentiating Eq. (jAll) with respect to r, we obtain 

a(t, r) 8'[r - r p (t)] + d r a(t, r)S[r- r p (t)] = a(t) S'[r - r p (t)]. 

Rearranging terms and using the rule of fully evaluating whenever possible, we find 

a(t,r)8'[r-r p (t)]=a(t)S'[r-r p (t)}-p(t)8[r-r p (t)}, p(t) = d r a(t,r p (t)) 
which is the analogous fully evaluated form. Upon integration, 

a(t,r) 8'[r ~ r p {t)\ dr = -p(t) = -d r a(t,r p (t)). 



d r a(t, r) 



r=r p (t) 



(A3) 



(A4) 



(A5) 



(A6) 



Since the first term on the right of Eq. (|A5[) disappears upon integration, why retain it? The answer is that we may 
multiply Eq. (| A5|) by another smooth (test) function, j(t, r). We can then proceed to fully evaluated form by reducing 
the smooth function j(t, r) a(t, r) on the left or use the same reduction on the first term on the right. In either case 
the result is 

7 (t, r) a(t, r) 6' [r - r p (t)] = j(t) a(t) S'[r ~ r p (t)} ~ a(t) d r j(t, r p {t)) S[r - r p {t)] - %t) d r a(t, r p {t)) S[r - r p (t)}. (A7) 

From this it is evident that we can partially evaluate a coefficient of 5' in a number of different ways. 

Martel [301 ) introduced the notation found in Eq. for gravitational master function source terms, with two- 

dimensional functions Gg m (t,r) and Fg m {t,r) multiplying 8 and 8', respectively. In examining the Zerilli-Moncrief 
master function, he left these coefficients partially evaluated. Sopuerta and Laguna [35j j started with the same notation 
for Gi m (t,r) and Fi m (t,r) in the case of the Cunningham-Price- Moncrief master function, and fully evaluated these 
coefficients at r = r p (t). A difficulty with the Gi m (t, r) and Ft m (t, r) notation is that there is no unique form of these 
functions if partially evaluated. Any solution of the RWZ wave equation will require a full evaluation of the source. 
The procedure should not matter but we prefer the clarity afforded by using the identities found in Eqs. (IA1[) and 
[5|l to write Eq. (jl.ll) in fully evaluated form from the outset 



S im (t,r) = G em {t) 8[r - r p (t)] + F tm {t) 8'[r - r p (t)], 



where 



Ge m (t) = G lm (t,r)- d r F tm (i, r) 



=r P (t) 



F lm {t)= F em (t,r) 



r P (t) 



(A8) 



(A9) 



Appendix B: Source terms for eccentric motion on Schwarzschild 

Here we give the unambiguous expressions for Gi m and Fg m for the even-parity Zerilli-Moncrief and odd-parity 
Cunningham-Price- Moncrief master functions fully evaluated at r — r p (t). We introduce new notation for constituent 
parts of Gim and Fg m based upon the projections of the stress-energy tensor defined in App.[C]and the fully evaluated 
time-dependent magnitudes of 8[r— r p (t)} given by Eqs. (I5.8[) and (|5.17l) . Note that we use Q and T to denote additional 
time-dependent factors that multiply the various stress-energy magnitudes. The indices on these Q and T factors are 
not tensor indices. 
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r 



1. Even parity 

In the even-parity case, we examine the terms first published by Martel (30j . but now fully evaluate them at 
r p (t). We find, 



Gim(t) = Q7 qVm + G? OL + Gi lira + G\ <?L + Q\ ?L 

Fi m (t) = T\ r <fg m + Tf 



~tt „tt 



(Bl) 



where 



(A + l)r p A2 



Gi r (t) = ^ , ^ A2 (A + f ) (Ar p + 6M) r p + 3M 



2 



r 



= " TYX^ A2 A (A + 1) Tp + 6AMr p + 15M 2 



(A + l)r p A2 



(B2) 



j „,2 3 



(A + 1)A P ' ^ W -(A + 1)A P ! 



with the q's given in Eq. Qb 



2. Odd parity 

In the odd-parity case, the fully evaluated source magnitudes are equivalent to those first published by Sopuerta 



and Laguna [35| and later with more detail by Field, Hesthaven, and Lau [44|. We find 



drf 

G tm (t) = QY P 7 em + G? -g* + G\p\ m , F im (t) = F t p\ m + T\ V \ m , (B3) 



where 



and the p's are given by Eq. (|5.17|) . 

Appendix C: Metric perturbation formalism in the Regge- Wheeler gauge 

Here we briefly summarize the definitions of metric perturbation (MP) amplitudes (on a common tensor spherical 
harmonic basis) for both even and odd parities. The field equations and Bianchi identities are given in terms of the 
MP amplitudes and spherical harmonic projected source terms. The specific gauge invariant master functions we use 
in our simulations are expressed in terms of the MP amplitudes and their associated master equations, potentials, 
and source terms are summarized. In what follows, lowercase Latin indices will run over (t, r), while uppercase Latin 
indices will run over (0,<p). This section draws heavily from Martel and Poisson |25j . The material here serves as a 
basis for discussing in Sec. IVlhow the MP can be numerically reconstructed from the master functions. 

1. Even parity 

Of the ten MP amplitudes, seven are in the even-parity sector. Using the decomposition of Martel and Poisson 
(25j . they are 

Pab (2H = 2^ h ab Y , PaB [ST) =2_^3 a Y B > PAB (^) = V ^ \K il AB Y + G Y AB j . (CI) 
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The tensor £Iab is the metric on the unit two-sphere, 

ds 2 = n AB dx A dx B = dd 2 + sin 2 9 dip 2 . (C2) 

The even-parity scalar (Y ), vector (Y^ m ), and tensor (Yi 1 ^ and VIabY ) spherical harmonics are defined in [25j. 
Note that Y^ S is the trace-free tensor spherical harmonic, which differs from what Regge and Wheeler used in their 
original work [19( • For the remainder of this section, we drop £ and m indices for the sake of brevity. 

In Schwarzschild coordinates, the amplitudes defined here are related to Regge and Wheeler's original quantities. 
In the u t,r sector," h u = fHo, h tr = Hi, and h rr = H%/f. For the off-diagonal elements, jt — ho and j r = hi. 
Finally, on the two-sphere Ghcrc = Grw, while inhere = i^RW — + l)G/2. We use the Regge- Wheeler gauge, where 
j a = G = 0. In this gauge and in Schwarzschild coordinates, the even-parity field equations are 

n2 „ 3r - 5M „ „ f _ , (A + 2) r + 2M A „ „ tt 

-d 2 K — — d r K + -d r h rr + i h rr + -^K = Q u , 



2 f 



2 f 

r-3M / A + l 
d t o r K H ^-—dtK d t h rr ^—h tr 



2 f 



-d 2 K- 



(r-M)f ^ 2/ / (A + l)r + 2M 

9 r i^ + —Othtr ~ -O r htt H ^3 ft-it - - 



/ 2 



9tft, rr — 9 r /i{ r + —dtK Tjhtr 

f r 2 f 

«», ^ai, MB- r ' M . (r - A/)/ 



r 2 / 



-a ( 2 /w + 2a t a r ^ tr - d 2 /i« - y9 t 2 if + /9 2 ^ + 2( - r r M) 



r-3Af 



(r - M )/ 



+ 



2(r — M) 



d r K 



(A + l)r 2 - 2(A + 2)Mr + 2M 
^T 2 



r2 



(A + l)r 2 - 2AMr - 2M 2 



jh tt - fh r 



Q r , 
Q\ 



(C3) 



which rely upon the following source terms 

Q ab {t, r) = 8vr J T ab Y* dfl, Q a {t, r) 

Q b (t, r) = 8nr 2 J T AB VL AB Y* dQ, Q J (t,r) = 32irr 
The conservation (Bianchi) identities are 



W+T) 



T aB Yg dn, 



(1-2)! 
(l + 2)\ 



(C4) 



t ab yx b dn. 



d t Q u + OrQ* 



2 f 



-Q tr - 



8 t Q tr + 8 r Q rr + ^-Q tt 



2r-5Af rr A + l 

r 2 f Q ~— ( 



= 0, 



(C5) 



A 



d t Q l + d r Q r + -Q r + Q"- — Q» = 



We use the gauge invariant Zerilli-Moncrief master function (see [13, HH, modifying the approach of [20]), which is 



* even (t,r) 



2r 



e(£ + i 

in Schwarzschild coordinates. It satisfies the wave equation 

d 2 d 2 



K+-{,f 2 h rr -rfd r K) 



dt 2 dr 2 



-V m 



= 5 e 



(C6) 



(C7) 



with source term 



Seven 



(A + l) A 



r 2 f (f 2 d r Q tt ~ d r Q rr ) + r(A - f)Q rr + r/ 2 Q b 

- ^ (A(A - l)r 2 + (4A - 9) Mr + 15A/ 2 ) Q u 



2/ 
A 



Q\ (C8) 
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and standard Zerilli potential 

V cvcn (r) 



r 2 A 2 



2A A + 1 



3M 
r 



18M 2 



. M 
A + — 

r 



(C9) 



2. Odd parity 

The remaining three MP amplitudes belong to the odd-parity sector, 

Pab = 0, p aS = £ ^a m ^i m , Pas (a*) - £ 



AB- 



(CIO) 



The vector (Xg) and tensor (X^g) spherical harmonics are those defined in [251 ] . Note that the tensor spherical 
harmonics differ from those used by Regge and Wheeler by a minus sign. For the remainder of this section, we again 
drop I and m indices. 

These MP amplitudes are related to Regge and Wheeler's quantities through h t — h , h r — hi, and h\ CIC = —h™. 
We use Regge- Wheeler gauge, in which h 2 — 0. In this gauge and in Schwarzschild coordinates, the odd-parity field 
equations are 



-d t d r h r + d 2 r h t - -d t h r 2{X + 1) 3 r - 4M h t = P\ 
r r° j 

2 2A / 

d?h r — d t d r h t H — d t h t H ^-h 

r r z 

1 2M 

--xd t h t + fd r h r H —h r = P. 

J H 



P r , 



(Cll) 



with source terms given by 



P a (t,r) 



IQitr 2 

The conservation (Bianchi) identity is 



T aB X* B dfl, 



P(t,r) = 167rr 4 ^4rT / T AB X* AB dn. 



2 2A 
dtP* + d r P r + -P r - ^rP = 0. 



(C12) 



(C13) 



In the odd- parity sector, we use the gauge-invariant Cunningham- Price- Moncrief master function [231 ] . which in 
Schwarzschild coordinates is 



$ odd (t,r) = - 



d r h t - d t h r h t 

r 



It satisfies the wave equation 



with source term 



dt 2 + dr 2 ° dd 



dd , 



rf 

S dd(t,r) = — 



-d t P r + fd r P l + ™P 
f r z 



and standard Regge- Wheeler potential 



V od d(r) = -o 



■cm)-** 

r 



(C14) 



(C15) 



(C16) 



(C17) 
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Appendix D: Asymptotic expansions for Jost functions at r* — > oo 

We examine here the asymptotic expansions that we use to set boundary conditions far from the black hole. The 
unit normalized solution to Eq. (|2.16[) is factored into the form 



where JeL„ is the "Jost function" [27|, which goes to 1 as 
function through R7 m = J, 



Irnn 



Eq. (|2.16[) and changing to r derivatives, we have 

"2M 



" mnr % (Dl) 

-oo. (We can similarly define the horizon side Jost 
^■ r * , which goes to 1 as r, -> —oo.) Plugging this into the source free version of 



/ 



d 2 7+ 



dr 2 



+ 2iu) n 



dr 



Vt , 

— t + — n 

j J £mn ~ u - 



(D2) 



From here we assume an asymptotic series solution of J^ mn of the form 



(D3) 



Note that contrary to a Taylor expansion which converges for fixed r with increasing j, this series converges for fixed 
j with increasing r. When a specific potential is chosen, the method of Frobenius can be used to find the coefficients 
dj . Plugging in the even-parity potential from Eq. (|C9[) a recurrence relation for the cij is 



2iX j aj — A 



A(i-l)i-12i<r(i-l)-2A(A+l) 



a,-i 



2(7 



A (3 - A) (j - 2) (j - 1) - (A 2 + 9ia) (j - 2) - 3A 5 



aj- 2 



3a' 



(3 - 4A) (j - 3) (j - 2) - 4A (j - 3) - 6A a,_ 3 - 18a 3 (j - 3) 2 a,_ 4 (D4) 



where u = Mui mn . For the odd- parity expansion, we plug in the potential in Eq. (|C17[) . The resulting recurrence 
relation is 



2ij cij 



-2a 



(J + 1) (J ~ 3) 



a,_2 



Z(£ + l)-j(j-l) 



aj-i- 



(D5) 



In order to use these recurrence relations, the first few terms ao, a\, . . . are needed. The recurrence relations actually 
provides them if one assumes that cij — for all negative j. 
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